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ABSTRACT 


In solar wind, magnetohydrodynamic (MHD) discontinuities are ubiquitous 
and often found to be at the origin of turbulence intermittency. They may also 
play a key role in the turbulence dissipation and heating of the solar wind. The 
tangential (TD) and rotational (RD) discontinuities are the two most important 
types of discontinuities. Recently, the connection between turbulence intermit- 
tency and proton thermodynamics has been being investigated observationally. 
Here we present numerical results from three-dimensional MHD simulation with 
pressure anisotropy and define new methods to identify and to distinguish TDs 
and RDs. Three statistical results obtained about the relative occurrence rates 
and heating effects are highlighted: (1) RDs tend to take up the majority of the 
discontinuities along with time; (2) the thermal states embedding TDs tend to 
be associated with extreme plasma parameters or instabilities, while RDs do not; 
(3) TDs have a higher average T as well as perpendicular temperature T';. The 
simulation shows that TDs and RDs evolve and contribute to solar wind heating 
differently. These results will inspire our understanding of the mechanisms that 


generate discontinuities and cause plasma heating. 


Subject headings: solar wind — magnetohydrodynamics — methods: numerical 
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1. Introduction 


The turbulent solar wind embodies discontinuities (e.g. m 
tangential discontinuity (TD) and rotational discontinuity (RD) are the two most important 
yet quite different types: in the deHoffmann-Teller frame, theoretically, TDs have no normal 
components of v and B on both sides, while RDs have normal components, must obey the 
Walén relation v = +B/,/fiop, and keep |v| and |B| continuous. Furthermore, TDs allow 
jumps of density and temperature, while these parameters have the same values on both 


sides of the RDs. As to the mechanisms of how discontinuities form, there exist two kinds of 
explanations. There are empirical evidences indicating that discontinuities are boundaries 
of magnetic flux tubes - ood T and there are suppositions that 
they form locally through nonlinear interactions and may be associated with small random 
currents mm m In Eu: 2D magnetohydrodynamic (MHD) 
simulation, most discontinuities appear to be TDs. However, in a 3D geometry, it remains 


unknown which type of discontinuity dominates. 


Dissipation of turbulence is considered an important contributor to the heating 
of the solar wind. Many recent studies concentrated on the role of intermittency and 
discontinuities in this process. |Bale et al. (2009) discovered strongly enhanced fluctuations 


Osman et all pou) 


2011)) reported that high PVI 


along the thresholds of plasma instabilities. 
(Partial Variance Increment) levels of various parameters correspond to intensive plasma 
heating and higher temperatures of electrons as well as ions. researched 
a large sample of data from measurements made by the Wind spacecraft and plotted the 
data distributions in the (8j, A) parameter plane (3) = p),/(B?/(2u0)), A = T, /Tj). Thus 


they could show that the distributions are roughly bounded by curves corresponding to the 


mirror and oblique fire-hose instabilities, that the regions near the instability thresholds 


m 


have higher averaged PVI, and that events with intense PVIs have A far from unity. 
(2013) analysed observed discontinuities with the PVI technique and found that 
the majority of them are RDs, but TDs have more obvious proton temperature increases. 


These empirical findings inspired us to investigate also the heating effects at the TDs and 


RDs obtained in our simulation. 


Numerical simulations have been employed before to understand plasma heating. 
(2008) assigned a path through the computational domain and then adopted 


the notion of PVI to analyse the 3D simulation data sampled along that path. Within 
their Hall MHD model they thus identified small-scale discontinuities being associated with 
intermittency. demonstrated by use of a 2.5D hybrid model that an 
ion temperature anisotropy can be created while the protons are heated by magnetic energy 
dissipation. conducted a full particle simulation which showed 
that, triggered by the Kelvin-Helmholtz instability with strong velocity shear, a turbulent 
cascade generates current sheets heating the plasma locally, and which yielded anisotropic 
particle distributions in that process. allowed for a broader range of 
Bj and the strength of the magnetic field fluctuations, thus obtaining results that basically 
are in accordance with those of Lu. However, neither were the data with 
high PVIs investigated, nor was the related heating of particles studied. All the mentioned 


work did also not distinguish between the various types of discontinuities in their simulation 


data. 


Motivated by all these aspects, we will here conduct a 3D numerical simulation with 
the aim to test new numerical methods to identify and analyse discontinuities, without 
assuming auxiliary paths along which data are sampled in the simulation domain. Based 
on this discontinuity identification, we present statistical results in order to investigate the 


proportion of TDs and RDs in all the discontinuities found, their parameter distribution in 
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the (8j, A) plane, and the temperature increases in TDs and RDs. Such counts of TDs and 
RDs and their distributions have, to our knowledge, not been reported previously. These 
new simulation results will help us to understand better particle heating at intermittent 


structures in the solar wind, and thus to resolve the turbulence dissipation problem. 


In Section 2 we will describe the numerical tools and the methods employed to identify 
and categorize discontinuities. In Section 3 we present a specific case of a TD and and RD, 
as well as statistical properties of all discontinuities found in the computation domain. In 


Section 4 we shall summarize our study, and further discuss relevant physical issues. 


2. Methods 
2.1. Numerical Model of MHD Turbulence 


In order to evaluate the role of temperature anisotropy in MHD turbulence, we adopt 


the model which employs the compressible ideal MHD equations and incorporates an 
anisotropic pressure tensor liens "mi : 


et 
2b iY. (v) =0, (1) 
8 I 
E pe (ovv t pil 4 (pj — »,)BB — (BB — B*1/2)/uo) Eu (2) 
PE +V x (-ux B) - 0 (3) 
Opi A n . Ópi 
we (pu) + 2pjB - (B - V)u = p (4) 
a "— 
E, + V- (pu) + 2p, V - u/34- 2(p) — p,)B - (B - V)u/3 = 0, (5) 


where p = (pij + 2p1)/3, B = B/ |B|. To describe instabilities correctly, it is common to use 


a heuristic term of pressure relaxation restricting the temperature anisotropy (Hesse & Birn 
1992; 1995 


denoted ópj/ót. In our numerical simulation we employ this 
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modified MHD model, which is not self-consistent as Vlasov models (e.g. 


b . This simplification of our model, on the other hand, will help to understand 
the role of fluid-like closures of the pressure tensor, which would be dissipative with the 
pressure-transport terms. Equations [4] and B] keep the quantities pı /(nB) and B?p, /n? 
adiabatically invariant ——! if the RHS of Equation [Jis zero. As Equation [3] 
is of the ideal MHD type, magnetic helicity is conserved. However, the model still lacks 


important kinetic physics of the solar wind, e.g. Landau damping and ion cyclotron 


resonances, which can be vital to turbulence dissipation. 


To solve the above system of equations, we utilize the BATSRUS codes (Powell et al. 
1999; [Tóth et a14/12012). The simulation is conducted in a three-dimensional cartesian 


coordinate system and encompasses an absolute volume of (62.8 Mm)? that is resolved 

in 256? grid points. The grid resolution (~ 250 km) is well above the ion skin depth 

~ 100 km), so Hall-dispersive physics is not included. We use the scheme proposed by 
in order to control numerical diffusion, and by applying such diffusion 
control strictly keep V - B — 0. This method guarantees proper dissipation and correct 
jump conditions at discontinuities. No explicit diffusive term is included in the numerical 
simulation code. Yet the magnetic and kinetic energy are subject to decay numerically. 
This decay can be attributed to the numerical scheme and grid resolution. We apply 
periodic boundary conditions to all the six surfaces of the simulation box. The initial 
conditions are set uniformly for p, T, and T,, and randomly for v and B, with the average 
Vo = 0 and a finite guide field Bo || é;. To simulate the solar wind at 1 AU, the field 


B, is chosen as 5 nT, while the Alfvén speed is set at 50 km s^! 


, and p = pj = p, with 
v 5p/(3p) = 50 km s~}, which corresponds to a proton number density of 5 cm? and 
an isotropic temperature of 10? K. For the turbulence part of the fields we take Fourier 


amplitudes obeying the broadband initial conditions described by [Matthaeus et al| (1996), 


with |óB(K)l^ = |dv(k)|? ox (1+ k/Knee) ? in the range 1077 m! < k < 8x 1077 m™. 
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The parameter ky nee is set to be 3 x 1077 m^! 


, and the spectral index q is set so that 
the slopes of the power spectral densities are both —5/3. We take ,/(v?) = 30 km s^! 


(accordingly 4/ (0B?) = 3 nT). The dimensionless cross-helicity o = 0.9. 


2.2. Methods of Selecting and Categorizing Discontinuities 


The analysis of discontinuities calls for reliable methods to select and categorize them 
in observational or numerical data. To trace abrupt spatial changes of the magnetic field, 
we use the total variance of increments (TVI) as an indicator and set a threshold for it. To 
calculate TVI, the total variance is first computed as 

JAB] =| M; (8Bs/00)?, (6) 
a,B=2,y,z 
where the partial derivative at grid point (i, j, k) about x is computed as 


B.(i +w, j, k) - B,(à — w, j,k) (7) 
2wóz l 

Here óx is the grid distance, B, denotes the corresponding component of magnetic field, 
and w is the width (in this work, we take tw — 3, i.e. within a cuboid of 7? grid points). 


For the y and z derivatives, similar differences along the corresponding directions are used. 


Then the TVI is the normalized total variance 
TVI = |AB| / (AB]), (8) 


where the angle brackets denote the average over the whole computational domain. This 


definition is a further development of the previous PVI defined by 2008 
(similar to the method adopted before by [Marsch & T m ), which was taken along a 
given path and hence directionally sensitive. Yet the above TVI includes all directions and 


thus is unbiased for all points possibly belonging to a discontinuity. The TVI utilises more 
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information from fully 3D data, and accordingly gives more physical insight. However, 
most solar wind measurements are 1D samples, so the method is inapplicable to such 
measurements. In the present work, those grid points with TVI > 3 are actually identifiable 


as discontinuity points and chosen for subsequent analysis. 


At each discontinuity point, we conduct the minimal variance analysis (MVA) 


in its neighbourhood, defined as a cuboid of the same given size 
for all points considered. Since Banı = B, is required, the direction of minimal variation 
can be regarded as the normal of the discontinuity EE E l In this work, 
we just consider a discontinuity locally as a small plane that contains the discontinuity 


point and cuts its neighbouring cuboid into two segments, so that averages of quantities on 


either side of the plane can be calculated in the segments obtained. 


To categorize the discontinuities, we use the criteria defined by (1973), which 


aim at judging two features: (1) whether B, = 0 holds, and (2) whether |B| remains 
continuous. Hence the parameters P, = |B,| / By, and P) = ô |B| / By are calculated (Bz is 
the larger of (|B|) on both sides). The points with P, < 0.2 and P, > 0.2 are categorized as 
TD, while the ones with P, > 0.4 and P5 < 0.2 as RD. All the other cases are categorized 
as ED(either TD or RD type, with P, < 0.4 and P, < 0.2) and ND(neither TD nor RD 
type, with P, > 0.2 and P» > 0.2). 


The aforementioned analysis only involves magnetic field data. To check and 
corroborate the results thus obtained, we also analyse the plasma velocity, density, and 
temperature data. In such case studies, it is trivial to check whether the density or 
temperature is continuous, but the velocity has to undergo the Walén test for RDs, i.e. one 
has to test whether v — vyr = +b = £B/,\/fop, which is usually done statistically in a 
scatter plot. The Walén test must be conducted in the deHoffmann- Teller frame. To find its 


velocity vyp, in the cuboid the sum $` jp) |(Vijk — Vo) X B;j,|^ is to be minimized. The 
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vo that makes this sum a minimum is then accepted as velocity of the deHoffmann-Teller 


frame Essen m 


3. Results 


The methods described above permit us to select the TDs and RDs from all cases with 
large TVI, and so we obtain corresponding data sets of discontinuities on which we can do 


individual and statistical research. 


To understand the evolution of the decaying turbulence in our MHD simulation, we plot 
in Figure[I]the squared current density (j?) averaged over the whole computational domain. 
This quantity relates to the curl of the magnetic field and describes its inhomogeneity 
and energy conversion (dissipation). Its evolution in time clearly shows the following 
phases: an initial drop, subsequent increase and final decay. This whole trend basically 
agrees with that found by in their previous simulation. To further 
illustrate that process, the z-component of current density j, is also plotted in a space-time 
display. The evolution implies that the initial drop of (j?) is due to a consumption of 
the magnetic energy in the compressive turbulent plasma motion, which leads to growing 
density inhomogeneity. This increase involves the formation of thin and stretched current 
sheets (see the numerous thin and sharp structures in (b3)). Due to the decay process, the 
inhomogeneities in (b4) fade away, yet a few current-sheet-like structures remain. The trace 


power spectral density of v, b, and power spectral density of n are also shown. A region 


with spectral index close to —5/3 seems to be identified. 


To check our methods and investigate the physical properties of the discontinuities, an 
individual TD and RD are picked and illustrated in Figure P| From the TD data we obtain 
|B„|/|B|z = 0.17 and ô| B|/|B|r, = 0.45. The TD has a normal ñ = (0.893, —0.355, —0.278), 
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almost aligned to the x-axis, and an HT velocity vyr = (4.13, —10.96, —41.99) km s !. 
In the HT frame, Panels (a) and (b) show that B and v across the discontinuity are 
quasi-parallel to the TD plane. The light grey cloud shows the sub-volume with high TVI 
and the plane is soaked therein. Panel (c) gives the TD's Walén test, where the points are 
rather scattered, but there still exists a correlation of v and b = B/ pop, especially in 
the y and z components (the correlation coefficients are Tss, ry, "zs = —0.06, 0.88, 0.82). 
For the RD we find |B,,|/|B|_ = 0.87, 9|B|/| B|r, = 0.11, n = (0.258, —0.139, 0.956) and 
vyr = (1.69,2.17, —55.56) km s !. Panels (d) and (e) have almost identical and slightly 
bent stream lines, thus illustrating the confirmation of the Walén relation. Panel (f) 
shows the correlation coefficients (Taz, yy, "zz = 0.87,0.97,0.86). The temperatures at 
the respective discontinuity points are also computed. The TD has T = 1.48 x 10° k, 
Tj = 1.51 x 10° K, and T, = 1.45 x 10° K, while the RD has T = 1.30 x 10? K, 

T| — 1.71 x 10? K, and T, — 1.09 x 10? K. The TD is by 13.596 hotter than the RD in T, 


and by 33.1% hotter in T. 


To emphasize furthermore the differences between TDs and RDs, statistical results are 
presented in Figure [8| where we plot the numbers of discontinuity points of each type as 
a function of time (left panel), as well as their percentages (right panel). At t — 0 there 
are only a few discontinuity points (53 TDs, 25 RDs, 65 EDs, 65 NDs; RDs are even not 
primary), but then the total number of discontinuity points increases with time, with RD 
becoming the dominant type. As the temporal evolution progresses, the number of TDs 
and their percentage first increase yet then decrease again slowly, behaving nearly in phase 


with that of the changing (j°), except during its initial drop phase. 


Moreover, in Figure [4] we plot for TDs and RDs their beta and anisotropy locations in 
the (8j, A) plane, where darker bins correspond to a higher number of discontinuity points, 


and with a uniform colour scale to facilitate comparisons for different times. For reference, 
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the threshold curves of the fire-hose and mirror instabilities (Hellinger et_al.||2006) are also 


plotted as orange and red curves, respectively. Apparently, the total numbers of TDs are 
less than those of RDs. Apart from that, they also distribute differently. At t — 60 s, the 
TD points tend to aggregate both in the centre region and near to the instability curves. 
In the decaying phase, the distribution shrinks toward the centre, and disappears finally. 
The RD points do not show this trend, with their majority still gathering around j = 0.5, 
A — 1.0, i.e. the initial values. At 60 s, some points lie beyond the instability lines but do 
not congregate there. For reference, the distributions of all grid points are supplied. They 


resemble those of the RDs. 


To investigate the heating effect of a discontinuity, the distributions of the temperatures 
found at the TD and RD points are also provided. Note that simulation of non-adiabatic 
heating is beyond the scope of this work as no real dissipation is involved. In Figure [5] we 
plot the distributions of T, T and T, at t = 60 s, with the values for TDs in red and RDs in 
blue. The TDs are slightly hotter in T and Tı. The TDs have T = 1.53 x 10? K (standard 
deviation c = 1.84 x 10* K), whereas the RDs have T = 1.30 x 10° K (c = 1.68 x 10* K). 
Also, TDs have T, = 1.42 x 10° K (e = 2.37 x 10* K), and RDs have T, = 1.19 x 10° K 
(c — 2.31 x 10* K). Since both types of discontinuity evolved in the same plasma with a 
uniform initial temperature, it is reasonable to conclude that TDs tend to become hotter 
than RDs, and to infer that TDs may intrinsically be more heated than RDs. The TDs' 
increased temperatures may be caused by plasma squeezing and adiabatic compression. In 


the present case study, the TD has a squeezing (quantified as V4,v4,) whose magnitude is 


5.41 times that of the RD. 
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4. Conclusion 


To conclude, in this letter we describe and apply methods to identify and analyse the 
ubiquitous discontinuities appearing in 3D simulation data, and can therefore present the 
following statistical results: (1) among the identified discontinuity points, RDs represent 
the majority, (2) TDs aggregate near the extremes of the fire-hose and mirror instability 


thresholds, while RDs do not, (3) TDs are hotter than RDs, both in their T and Tj. 


However, this work is also limited in some aspects. Though our simulation of decaying 
turbulence in a closed box can reveal its different stages of evolution in time, the results 
obtained are difficult to compare with turbulence observations made in the solar wind, 
where the convected plasma volume is open and can always be filled with fresh waves 
continuously supplied from a source region. To alleviate this problem (i.e. the difference 
between the reality and our case study and its temperature distributions), we selected the 
results obtained at the time when (j?) is maximal, which may represent a state being close 


to developed turbulence. 


Moreover, MHD though with thermal anisotropy lacks realistic descriptions of 
microscopic (dissipative) solar wind processes. However, the used MHD model allows 
us to reduce significantly the computational cost and to investigate discontinuities in 
three-dimensional space, a virtue of our approach which appears crucial for the analysis. 
In the identification of discontinuities, we take w — 3 and the TVI threshold as 3, values 
which are reasonably chosen but seem somehow arbitrary. Therefore, we also checked the 
outcome with a different w (ranging from 2 to 5) and TVI threshold (set lower at 2.64) 
for the statistical study, but found that the results did not change essentially, in terms of 


identification, normal direction, proportion, and distribution of all discontinuities. 


The physical mechanisms generating TDs and RDs should be further investigated. 


From another aspect, hybrid or full particle simulations should be considered, as they 
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enable the description of kinetic processes in the solar wind as well. Besides, even in 3D 
MHD simulations, the methods to study large TVI events could be further improved, e.g., 
possibly by identifying the discontinuities with model structures having a 3D geometric 
configuration instead of by simply counting their associated points. Such approach may 
improve the identification of the discontinuities and give better statistics concerning their 


occurrence rates as well as the percentages of their local temperature increase. 


This work is supported by NSFC grants under contracts 41231069, 41174148, 
41222032, 41274132, 41474147, 41031066 and 41304133, and was carried out using the 
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Fig. 1.— Evolution stages of the decaying turbulence in the MHD simulation. (a): Temporal 
profile of (7?) in arbitrary unit. (b): Distribution of the z-component of the current density 
jz in the plane z = 31.4 Mm at the given times; 7, is given in arbitrary unit in agreement 
with (a). (c): Trace power spectral density of v (red) and b (blue), as well as PSD of n 
(green) at the same times. The horizontal axis represents k, = \/k? + K2, while the vertical 
axis represents v or b(black, left) and n(green, right). Spectral indices —5/3(black) and 


—A(cyan) are also plotted in dashed lines for reference. 
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Fig. 2.— A TD (top panels) and RD (bottom panels) at t = 60 s. (a) and (b): Schematic of 
the TD, where the transparent blue plane shows the plane of discontinuity, the red and yellow 
arrows denote respectively B and v, the light grey clouds plot the TVI with thicker cloud 
representing larger values, and the red, blue and green arrows at the box corner indicate 
the x, y and z directions, respectively. The background guide field is along the z-direction. 
(c): Scatter plot of the components of v’ = v — vyr and b = B/ nop of the TD. (d) and 
(e): Schematic of the RD displayed with stream lines; cyan arrows show the direction of the 


velocity in the HT-frame. (f): Scatter plot for the RD. 
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Fig. 3.— Counts (left) of discontinuity points and their proportions (right) for each type. 


TDs, NDs, RDs and EDs are represented respectively in red, green, cyan and purple. 
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Fig. 4.— Temporal evolution of the distributions of TD points (top), RD points (centre), and 
all grid points in the whole computational region (bottom) in the (5j, A = Tı /Tj) plane, 
arranged in columns for different times. For reference, the red and orange lines give the 
thresholds for the mirror and fire-hose instabilities, respectively. The colour of a bin denotes 


the number of grid points therein. 
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Fig. 5.— Distributions of the temperatures associated with the TD and RD points. The 


left, middle and right sub-plot shows the distribution of T, Tj and T! , respectively. 


